# =============================================================================
# Title: Path Analysis Model 1 - SBL, GC, JC, and Work Procrastination
# Description: Moderated parallel mediation analysis examining how supervisor 
#              bottom-line mentality affects work procrastination through 
#              goal clarity and job calling, with red tape as moderator
# =============================================================================

# Load required packages
library(lavaan)  # for path analysis
library(readxl)  # for reading Excel files
library(psych)   # for descriptive statistics

# Function for Monte Carlo confidence intervals
monte_carlo_ci <- function(a, b, se_a, se_b, rep = 20000, conf = 0.95) {
  set.seed(123)  # for reproducibility
  a_dist <- rnorm(rep, mean = a, sd = se_a)
  b_dist <- rnorm(rep, mean = b, sd = se_b)
  indirect <- a_dist * b_dist
  
  ci_lower <- quantile(indirect, (1 - conf) / 2)
  ci_upper <- quantile(indirect, 1 - (1 - conf) / 2)
  mean_indirect <- mean(indirect)
  se_indirect <- sd(indirect)
  
  return(c(estimate = mean_indirect, 
           se = se_indirect,
           ci.lower = ci_lower, 
           ci.upper = ci_upper))
}

# Function to calculate R² significance and SE
calculate_r2_statistics <- function(r2, n, k) {
  f_stat <- (r2 / k) / ((1 - r2) / (n - k - 1))
  p_value <- 1 - pf(f_stat, k, n - k - 1)
  se_r2 <- sqrt((4 * r2 * (1 - r2)^2 * (n - k - 1)^2) / ((n^2 - 1) * (n + 3)))
  return(list(f_stat = f_stat, p_value = p_value, se = se_r2))
}

# -----------------------------------------------------------------------------
# Data Preparation
# -----------------------------------------------------------------------------

# Read the Excel file
data <- read_excel("C:/Users/wangm/Nutstore/1/AI & Big Data Group/领导底线心智BLM-工作拖延/Supplemental materials/CC/Supervisor Bottom-Line Mentality_Recoded.xlsx")


# Center predictors for interaction analysis
data$SBL_c <- scale(data$SBL, scale=FALSE)
data$RT_c <- scale(data$RT, scale=FALSE)
data$SBLxRT <- data$SBL_c * data$RT_c

# -----------------------------------------------------------------------------
# Descriptive Statistics and Correlations
# -----------------------------------------------------------------------------

# Calculate descriptive statistics
desc_stats <- describe(data[c("SBL", "RT", "GC", "JC", "WP", "CS", "HS",
                              "Gender", "Edu", "Age", "SStenure", "Rank")])

# Calculate correlation matrix
cor_matrix <- cor(data[c("SBL", "RT", "GC", "JC", "WP", "CS", "HS",
                         "Gender", "Edu", "Age", "SStenure", "Rank")])

# -----------------------------------------------------------------------------
# Path Analysis Model Specification
# -----------------------------------------------------------------------------

# Define the path model
model <- '
  # Direct effects on mediators
  GC ~ a1*SBL_c + w1*RT_c + x1*SBLxRT + 
       c1*Gender + c2*Edu + c3*Age + c4*SStenure + c5*Rank +
       c16*CS + c17*HS  # Added CS and HS controls

  JC ~ a2*SBL_c + w2*RT_c + x2*SBLxRT + 
       c6*Gender + c7*Edu + c8*Age + c9*SStenure + c10*Rank +
       c18*CS + c19*HS  # Added CS and HS controls

  # Direct effects on outcome
  WP ~ c*SBL_c + b1*GC + b2*JC + x3*SBLxRT +
       c11*Gender + c12*Edu + c13*Age + c14*SStenure + c15*Rank +
       c20*CS + c21*HS  # Added CS and HS controls

  # Indirect effects
  ie1 := a1*b1    # SBL -> GC -> WP
  ie2 := a2*b2    # SBL -> JC -> WP
  total := c + (a1*b1) + (a2*b2)  # Total effect

  # Add residual covariances between mediators
  GC ~~ JC
'

# Function to format p-values
format_p <- function(p) {
  if (p < .001) return("< .001")
  if (p < .01) return(sprintf("%.3f", p))
  return(sprintf("%.3f", p))
}

# Function to add significance stars
add_stars <- function(p) {
  if (p < .001) return("***")
  if (p < .01) return("**")
  if (p < .05) return("*")
  return("")
}

# Estimate the model
fit <- sem(model, 
           data = data,
           estimator = "MLR",    # Robust maximum likelihood estimation
           se = "robust",        # Robust standard errors
           cluster = "organizationcode",  # Clustering by organization
           missing = "fiml")     # Handle missing data

results <- parameterEstimates(fit)
std_results <- standardizedSolution(fit)

# Create output file
sink("path_analysis_model1_results.txt")

# Print header
cat("\n============================================================================\n")
cat("                         Path Analysis Results - Model 1                      \n")
cat("============================================================================\n\n")

# Function to print model results
print_model_results <- function(dv, results, std_results) {
  cat(sprintf("\nDependent Variable: %s\n", dv))
  cat("----------------------------------------------------------------------------\n")
  
  if(dv == "WP") {
    cat(sprintf("%-15s %-10s %-10s %-10s %-10s %-10s\n", 
                "Variables", "Coef", "SE", "Std.Coef", "t", "p"))
  } else {
    cat(sprintf("%-15s %-10s %-10s %-10s %-10s %-10s\n", 
                "Variables", "γ", "SE", "Std.γ", "t", "p"))
  }
  cat("----------------------------------------------------------------------------\n")
  
  dv_results <- results[results$lhs == dv & results$op == "~",]
  std_dv_results <- std_results[std_results$lhs == dv & std_results$op == "~",]
  
  for(i in 1:nrow(dv_results)) {
    stars <- add_stars(dv_results$pvalue[i])
    beta <- sprintf("%.3f", std_dv_results$est.std[i])
    t_value <- dv_results$est[i] / dv_results$se[i]
    
    cat(sprintf("%-15s %-10.3f %-10.3f %-10s %-10.3f %-10s %s\n",
                dv_results$rhs[i],
                dv_results$est[i],
                dv_results$se[i],
                beta,
                t_value,
                format_p(dv_results$pvalue[i]),
                stars))
  }
  
  n <- nrow(data)
  k <- nrow(dv_results)
  r2 <- inspect(fit, "r2")[dv]
  r2_stats <- calculate_r2_statistics(r2, n, k)
  
  cat("----------------------------------------------------------------------------\n")
  cat(sprintf("R² = %.3f (SE = %.3f), F(%d, %d) = %.3f, p %s %s\n", 
              r2, r2_stats$se, k, n-k-1, r2_stats$f_stat, 
              format_p(r2_stats$p_value),
              add_stars(r2_stats$p_value)))
  cat("----------------------------------------------------------------------------\n")
}

# Print results for each dependent variable
print_model_results("GC", results, std_results)
print_model_results("JC", results, std_results)
print_model_results("WP", results, std_results)

# Calculate and print indirect effects
a1 <- results$est[results$lhs == "GC" & results$rhs == "SBL_c"]
se_a1 <- results$se[results$lhs == "GC" & results$rhs == "SBL_c"]
b1 <- results$est[results$lhs == "WP" & results$rhs == "GC"]
se_b1 <- results$se[results$lhs == "WP" & results$rhs == "GC"]

a2 <- results$est[results$lhs == "JC" & results$rhs == "SBL_c"]
se_a2 <- results$se[results$lhs == "JC" & results$rhs == "SBL_c"]
b2 <- results$est[results$lhs == "WP" & results$rhs == "JC"]
se_b2 <- results$se[results$lhs == "WP" & results$rhs == "JC"]

# Calculate Monte Carlo CIs
mc_ie1 <- monte_carlo_ci(a1, b1, se_a1, se_b1)
mc_ie2 <- monte_carlo_ci(a2, b2, se_a2, se_b2)

# Print indirect and total effects
cat("\nIndirect and Total Effects (Monte Carlo Method, 20,000 replications)\n")
cat("----------------------------------------------------------------------------\n")
cat(sprintf("%-20s %-10s %-10s %-10s %-10s\n",
            "Effect", "Estimate", "SE", "95% CI LL", "95% CI UL"))
cat("----------------------------------------------------------------------------\n")

cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
            "IE1 (SBL→GC→WP)",
            mc_ie1[1], mc_ie1[2], mc_ie1[3], mc_ie1[4]))

cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
            "IE2 (SBL→JC→WP)",
            mc_ie2[1], mc_ie2[2], mc_ie2[3], mc_ie2[4]))

# Print total effect
total_effect <- results$est[results$label == "total"]
total_se <- results$se[results$label == "total"]
cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
            "Total Effect",
            total_effect,
            total_se,
            total_effect - 1.96*total_se,
            total_effect + 1.96*total_se))

sink()

# Calculate conditional indirect effects
RT_sd <- sd(data$RT)

# Function for conditional indirect effects
conditional_indirect_effects <- function(results, sd_RT) {
  # Extract all necessary coefficients and standard errors
  a1 <- results$est[results$lhs == "GC" & results$rhs == "SBL_c"]
  se_a1 <- results$se[results$lhs == "GC" & results$rhs == "SBL_c"]
  x1 <- results$est[results$lhs == "GC" & results$rhs == "SBLxRT"]
  se_x1 <- results$se[results$lhs == "GC" & results$rhs == "SBLxRT"]
  b1 <- results$est[results$lhs == "WP" & results$rhs == "GC"]
  se_b1 <- results$se[results$lhs == "WP" & results$rhs == "GC"]
  
  a2 <- results$est[results$lhs == "JC" & results$rhs == "SBL_c"]
  se_a2 <- results$se[results$lhs == "JC" & results$rhs == "SBL_c"]
  x2 <- results$est[results$lhs == "JC" & results$rhs == "SBLxRT"]
  se_x2 <- results$se[results$lhs == "JC" & results$rhs == "SBLxRT"]
  b2 <- results$est[results$lhs == "WP" & results$rhs == "JC"]
  se_b2 <- results$se[results$lhs == "WP" & results$rhs == "JC"]
  
  sink("conditional_indirect_effects_model1.txt")
  
  cat("\n============================================================================\n")
  cat("                     Conditional Indirect Effects Analysis                    \n")
  cat("============================================================================\n\n")
  
  RT_levels <- c(-sd_RT, sd_RT)
  RT_labels <- c("Low RT", "High RT")
  
  for(i in 1:length(RT_levels)) {
    RT <- RT_levels[i]
    cat(sprintf("\nAt %s (RT = %.2f SD):\n", RT_labels[i], RT))
    cat("----------------------------------------------------------------------------\n")
    
    # Calculate conditional effects
    cond_a1 <- a1 + x1 * RT
    cond_a2 <- a2 + x2 * RT
    
    # Calculate standard errors
    se_cond_a1 <- sqrt(se_a1^2 + RT^2*se_x1^2)
    se_cond_a2 <- sqrt(se_a2^2 + RT^2*se_x2^2)
    
    # Calculate conditional indirect effects
    ie1 <- monte_carlo_ci(cond_a1, b1, se_cond_a1, se_b1)
    ie2 <- monte_carlo_ci(cond_a2, b2, se_cond_a2, se_b2)
    
    cat(sprintf("First-stage moderation for GC path (a1 + x1*RT): %.3f\n", cond_a1))
    cat(sprintf("First-stage moderation for JC path (a2 + x2*RT): %.3f\n", cond_a2))
    cat("\nConditional indirect effects:\n")
    cat("----------------------------------------------------------------------------\n")
    cat(sprintf("%-20s %-10s %-10s %-10s %-10s\n",
                "Path", "Effect", "SE", "95% CI LL", "95% CI UL"))
    cat("----------------------------------------------------------------------------\n")
    
    cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
                "SBL→GC→WP",
                ie1[1], ie1[2], ie1[3], ie1[4]))
    
    cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
                "SBL→JC→WP",
                ie2[1], ie2[2], ie2[3], ie2[4]))
  }
  
  # Index of moderated mediation
  imm1 <- monte_carlo_ci(x1, b1, se_x1, se_b1)
  imm2 <- monte_carlo_ci(x2, b2, se_x2, se_b2)
  
  cat("\n============================================================================\n")
  cat("                      Index of Moderated Mediation                           \n")
  cat("============================================================================\n")
  cat("----------------------------------------------------------------------------\n")
  cat(sprintf("%-20s %-10s %-10s %-10s %-10s\n",
              "Path", "Index", "SE", "95% CI LL", "95% CI UL"))
  cat("----------------------------------------------------------------------------\n")
  
  cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
              "Through GC",
              imm1[1], imm1[2], imm1[3], imm1[4]))
  
  cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
              "Through JC",
              imm2[1], imm2[2], imm2[3], imm2[4]))
  
  sink()
}

# Run conditional indirect effects analysis
conditional_indirect_effects(results, RT_sd)



# 创建格式化的结果表格文件
sink("formatted_path_analysis_table.txt")

cat("============================================================================\n")
cat("                         Path Analysis Results - Model 1                     \n")
cat("============================================================================\n\n")

# 打印GC模型结果
cat("\nDependent Variable: GC\n")
cat("----------------------------------------------------------------------------\n")
cat(sprintf("%-15s %-10s %-10s %-10s %-10s %-10s\n", 
            "Variables", "γ", "SE", "Std.γ", "t", "p"))
cat("----------------------------------------------------------------------------\n")

gc_results <- results[results$lhs == "GC" & results$op == "~",]
std_gc_results <- std_results[std_results$lhs == "GC" & std_results$op == "~",]

for(i in 1:nrow(gc_results)) {
  stars <- add_stars(gc_results$pvalue[i])
  beta <- sprintf("%.3f", std_gc_results$est.std[i])
  t_value <- gc_results$est[i] / gc_results$se[i]
  
  cat(sprintf("%-15s %-10.3f %-10.3f %-10s %-10.3f %-10s %s\n",
              gc_results$rhs[i],
              gc_results$est[i],
              gc_results$se[i],
              beta,
              t_value,
              format_p(gc_results$pvalue[i]),
              stars))
}

n <- nrow(data)
k <- nrow(gc_results)
r2_gc <- inspect(fit, "r2")["GC"]
r2_stats_gc <- calculate_r2_statistics(r2_gc, n, k)

cat("----------------------------------------------------------------------------\n")
cat(sprintf("R² = %.3f (SE = %.3f), F(%d, %d) = %.3f, p %s %s\n", 
            r2_gc, r2_stats_gc$se, k, n-k-1, r2_stats_gc$f_stat, 
            format_p(r2_stats_gc$p_value),
            add_stars(r2_stats_gc$p_value)))
cat("----------------------------------------------------------------------------\n")

# 打印JC模型结果
cat("\nDependent Variable: JC\n")
cat("----------------------------------------------------------------------------\n")
cat(sprintf("%-15s %-10s %-10s %-10s %-10s %-10s\n", 
            "Variables", "γ", "SE", "Std.γ", "t", "p"))
cat("----------------------------------------------------------------------------\n")

jc_results <- results[results$lhs == "JC" & results$op == "~",]
std_jc_results <- std_results[std_results$lhs == "JC" & std_results$op == "~",]

for(i in 1:nrow(jc_results)) {
  stars <- add_stars(jc_results$pvalue[i])
  beta <- sprintf("%.3f", std_jc_results$est.std[i])
  t_value <- jc_results$est[i] / jc_results$se[i]
  
  cat(sprintf("%-15s %-10.3f %-10.3f %-10s %-10.3f %-10s %s\n",
              jc_results$rhs[i],
              jc_results$est[i],
              jc_results$se[i],
              beta,
              t_value,
              format_p(jc_results$pvalue[i]),
              stars))
}

k_jc <- nrow(jc_results)
r2_jc <- inspect(fit, "r2")["JC"]
r2_stats_jc <- calculate_r2_statistics(r2_jc, n, k_jc)

cat("----------------------------------------------------------------------------\n")
cat(sprintf("R² = %.3f (SE = %.3f), F(%d, %d) = %.3f, p %s %s\n", 
            r2_jc, r2_stats_jc$se, k_jc, n-k_jc-1, r2_stats_jc$f_stat, 
            format_p(r2_stats_jc$p_value),
            add_stars(r2_stats_jc$p_value)))
cat("----------------------------------------------------------------------------\n")

# 打印WP模型结果
cat("\nDependent Variable: WP\n")
cat("----------------------------------------------------------------------------\n")
cat(sprintf("%-15s %-10s %-10s %-10s %-10s %-10s\n", 
            "Variables", "Coef", "SE", "Std.Coef", "t", "p"))
cat("----------------------------------------------------------------------------\n")

wp_results <- results[results$lhs == "WP" & results$op == "~",]
std_wp_results <- std_results[std_results$lhs == "WP" & std_results$op == "~",]

for(i in 1:nrow(wp_results)) {
  stars <- add_stars(wp_results$pvalue[i])
  beta <- sprintf("%.3f", std_wp_results$est.std[i])
  t_value <- wp_results$est[i] / wp_results$se[i]
  
  cat(sprintf("%-15s %-10.3f %-10.3f %-10s %-10.3f %-10s %s\n",
              wp_results$rhs[i],
              wp_results$est[i],
              wp_results$se[i],
              beta,
              t_value,
              format_p(wp_results$pvalue[i]),
              stars))
}

k_wp <- nrow(wp_results)
r2_wp <- inspect(fit, "r2")["WP"]
r2_stats_wp <- calculate_r2_statistics(r2_wp, n, k_wp)

cat("----------------------------------------------------------------------------\n")
cat(sprintf("R² = %.3f (SE = %.3f), F(%d, %d) = %.3f, p %s %s\n", 
            r2_wp, r2_stats_wp$se, k_wp, n-k_wp-1, r2_stats_wp$f_stat, 
            format_p(r2_stats_wp$p_value),
            add_stars(r2_stats_wp$p_value)))
cat("----------------------------------------------------------------------------\n")

# 打印间接效应
cat("\nIndirect and Total Effects (Monte Carlo Method, 20,000 replications)\n")
cat("----------------------------------------------------------------------------\n")
cat(sprintf("%-20s %-10s %-10s %-10s %-10s\n",
            "Effect", "Estimate", "SE", "95% CI LL", "95% CI UL"))
cat("----------------------------------------------------------------------------\n")

cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
            "IE1 (SBL→GC→WP)",
            mc_ie1[1], mc_ie1[2], mc_ie1[3], mc_ie1[4]))

cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
            "IE2 (SBL→JC→WP)",
            mc_ie2[1], mc_ie2[2], mc_ie2[3], mc_ie2[4]))

# 打印总效应
total_effect <- results$est[results$label == "total"]
total_se <- results$se[results$label == "total"]
cat(sprintf("%-20s %-10.3f %-10.3f %-10.3f %-10.3f\n",
            "Total Effect",
            total_effect,
            total_se,
            total_effect - 1.96*total_se,
            total_effect + 1.96*total_se))

sink()

cat("\n格式化表格已保存为: formatted_path_analysis_table.txt\n")